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This is a longer version of our article W, containing more detailed 
explanations and providing pedagogical introductions to the methods we 
use. 

We consider a product of an arbitrary number of independent rect- 
angular Gaussian random matrices. We derive the mean densities of its 
eigenvalues and singular values in the thermodynamic limit, eventually ver- 
ified numerically. These densities are encoded in the form of the so-called 
M-transforms, for which polynomial equations are found. We exploit the 
methods of planar diagrammatics, enhanced to the non-Hermitian case, 
and free random variables, respectively; both are described in the appen- 
dices. As particular results of these two main equations, we find the singular 
behavior of the spectral densities near zero. Moreover, we propose a finite- 
size form of the spectral density of the product close to the border of its 
eigenvalues' domain. Also, led by the striking similarity between the two 
main equations, we put forward a conjecture about a simple relationship 
between the eigenvalues and singular values of any non-Hermitian random 
matrix whose spectrum exhibits rotational symmetry around zero. 

PACS numbers: 02.50.Cw (Probability theory), 02.70.Uu (Applications of 
Monte Carlo methods), 05.40.Ca (Noise) 

1. Introduction and the Main Results 

1.1. Introduction 
1.1.1. Products of Random Matrices 

The problem of multiplication of random matrices received a consider- 
able measure of attention from the random matrix theory (RMT) commu- 
nity recently I2l[3llll[5l[6l[71[8l[9l[l0l[lll[l2l[l3]. In this matter, one is 
interested in studying the statistical properties (in particular, the average 
distribution of the eigenvalues) of the product 

P = AiA2...Ai (1) 

of L > 1 random matrices Ai, I = 1,2, . . . , L. In the basic version of the 
problem, all these matrices are assumed to be statistically independent, to 
which condition we too will adhere. 

In order for these matrices to be multiplicable, the dimensions of each 
have most generally to be of the form A'^; x Ni^i, where Ni, N2, ■ ■ ■ , Nl, Nl^i 
are integers. Then, P has sizes A^i x N^^i; if this matrix is to have eigen- 
values, it must obviously be square, i.e., N^-^^i = Ni. 

* Presented by Giacomo Livan at the 23rd Marian Smoluchowski Symposium on Sta- 
tistical Physics "Random Matrices, Statistical Physics and Information Theory," 
September 26-30, 2010, Krakow, Poland. 

^ Igiacomo . livanSpv. infn. it (corresponding author) 
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Moreover, in RMT one typically considers the "thermodynamic limit" in 
which the matrices are infinitely large, but their dimensions have the same 
order of magnitude, i.e., their ratios stay finite; here, we decide to define 
these ratios w.r.t. A''^^!, 

Ni 

Ni oo, with Ri = — = finite, for / = 1, 2, . . . , L, L + 1. (2) 

(Obviously, Rl+i = 1; if we additionally take NiJ^i = Ni, then there also is 
Ri = 1.) In this limit, all finite-size corrections are lost, such as universal 
oscillations of the spectrum at its edges. 



1.1.2. Products of Square Girko— Ginibre Random Matrices 

To date, only various ensembles of square {i.e., Ni = N2 = ■ ■ ■ = N^^i = N) 
random matrices have been investigated in the context of the multiplication 
task. For instance, the authors of [2] considered each A; to consists of 
complex entries whose real and imaginary parts are all IID Gaussian ran- 
dom variables with zero mean and variance af /2N; in other words, such 
a random matrix model, called the "Girko-Ginibre ensemble" [15\ [T6| [T7] , 
has the probability measure, 

d^(AOoce ^ ^ ' ^DA,, (3) 

where the flat measure DAj = n^b=i d(Re[A/]af,)d(Im[A/]afe), and the nor- 
malization constant has been omitted. Using the method of planar diagrams 
and Dyson-Schwinger's equations, which is the same technique we will em- 
ploy here in section [2] and pedagogically describe in appendix [XJ one may 
derive ( ^A.2.2p that on average, the eigenvalues of a Girko-Ginibre matrix 
are scattered within a centered circle of radius ai with the uniform density, 

. -^ f for lAI <ai. 

Employing non~Hermitian planar diagrammatics again, yet in a more 
demanding situation, the authors of [H] (see their equation (5)) discovered 
that the eigenvalues of the product P ([1]) of such Girko-Ginibre matrices 
also fill a centered circle, of radius 

(T = 0-1 (T2 •• -CTL, (5) 

with the average density given by a surprisingly simple expression, 

/5P(A,A) = / i;^l^r'^'~*^ for |A|<cT, (6) 
[ 0, for |A| >a. ^ ^ 
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Remarkably, this formula remains valid even if the constituent matrices are 
not identically distributed, just with the assumptions of independence and 
Gaussianity retained, i.e., they may come from different Gaussian ensem- 
bles, such as GUE, GOE, or the so-called "Gaussian elliptic ensembles." 
Moreover, it has been conjectured [14j that ^ holds for an even wider class 
of matrices, such as ones having independent entries fulfilling the Pastur- 
Lindeberg condition — in this sense, the result ^ is universal. One unex- 
pected implication of this universality is that a product of random matri- 
ces whose spectra do not necessarily display rotational symmetry has the 
eigenvalue distribution on the complex plane which does possess rotational 
symmetry (i.e., the average density depends only on |A|). 

1.1.3. Non— Hermitian Random Matrix Models 

Exploration of products of matrix models brings us almost inevitably 
into the realm of non-Hermitian random matrices. The most pronounced 
difference between them and their Hermitian cousins is that their eigenvalues 
are generically complex, while Hermitian spectra must be real. In the ther- 
modynamic limit, the eigenvalues of non-Hermitian ensembles cover two- 
dimensional domains on the complex plane, contrary to one-dimensional 
cuts in the Hermitian case. 

Much of the RMT machinery — such as the saddle-point method, or- 
thogonal polynomials, the Efetov's super symmetric technique, or the dia- 
grammatic expansion and the free random variables calculus (see for exam- 
ple [l8l[l9l[20l[2ll|22l|23llH[25l[M[2Zl[2g reviews; these last 
two techniques are used in this paper) — has been developed in the Her- 
mitian context, and tailored for real spectra. The necessity of dealing with 
complex spectra demands enhancements of these methods. In appendix \X\ 
we present a self-contained crash course on how the planar diagrammatics 
can be applied to find average spectra of non-Hermitian random ensem- 
bles taaiMiiM]. 

Leaving thus all the details for later, let us just mention now that as the 
information about the average spectrum of a Hermitian matrix model H is 
encoded in the "spectral density" ph(A) (|A.ip — which is a function of a real 
variable, or equivalently in the "Green's function" Gh(-z) (|A.3p — which is 
a holomorphic function everywhere except the cuts where the eigenvalues 
lie, or equivalently in the so-called "M-transform," M-ii{z) = zGii{z) — 1 
(IXTI) . so for a non-Hermitian matrix X, one exploits an analogous set of 
concepts, namely, the spectral density /Ox(A, A) (jA.lSp — which is now a 
function of a complex argument, denoted here by A and A as independent 
variables (this is the object appearing in ^ and ([6])), the "non-holomorphic 
Green's function," Gx(-2,^) ()A.20p . as well as the "non-holomorphic M~ 
transform," M-x_{z,'z) = zG^{z,'z) — 1 ()A.24p . Presentation in appendix 1X1 
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of the diagrammatic method augmented to handle non~Hermitian products 
is one of the motivations of this article. 

1.1.4. Singular Values and Free Random Variables 

Besides eigenvalues, another important characteristics of a matrix model 
P is given by its "singular values," which are defined as the square roots of 
the (real) eigenvalues of the Hermitian matrix 



We will investigate the singular values in section [3l using and thereby 
promoting an approach called the "free random variables" (FRV) calcu- 
lus. We sketch its fundamental implications in appendix [Bl The idea 
will be to rewrite Q, through cyclic permutations of the constituent ma- 
trices, as a product of certain Hermitian ensembles. Now, for a Hermi- 
tian product of independent (more precisely, "free" ()B.6|) — this is a gen- 
eralized notion of statistical independence, suited for matrix probability 
theory) Hermitian matrices — FRV provides a "multiplication algorithm" : 
Given two Hermitian random matrices. Hi and H2 whose product is still 
Hermitian, and which are mutually free, one begins from finding their 
"A'^-transforms" ()B.19|) . defined as the functional inverses of the respec- 
tive M-transforms, A^Hi.2(-^Hi.2(-2)) = AfHi,2(-^Hi,2(-^)) = ^- Then, a ba- 
sic FRV theorem claims that the A^-transform of the product H1H2 is 
in fact the product of the two A/'-transforms, up to a simple prefactor, 
N-R^Uiiz) = {z/{z + l))NYi^{z)N-ti2{z) dElQl). Inverting functionally the 
result, one obtains the M-transform, and thus the entire spectral informa- 
tion about H1H2. 



1.2.1. The Aim of the Paper 

In this paper, our objective is to generalize the work of [U] by calculat- 
ing both the eigenvalues and singular values of a product ([1]) of rectangular 
matrices with complex entries being IID Gaussian random variables. More 
precisely, for our L random matrices A; we will now allow arbitrary rectan- 
gularity ratios Ri ^ , while the assumption of the real and imaginary parts 
of all the matrix elements of each A; being IID Gaussian random numbers 
will be realized by taking the following probability measures. 



Q = Ptp. 



(7) 



1.2. The Main Results 





It is almost identical to ([3]), with the exception of the scaling of the variance, 
whose inverse has generically to be chosen of the same order of magnitude 
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as the dimensions of the matrices in question — in order to yield in the 
thermodynamic hmit a meaningful density of the eigenvalues — but oth- 
erwise it can be arbitrary; here, we decide for the factor of We 
will assume A'^^^+i = -^i (^-e., R\ = 1) when considering the eigenvalues of 
P Q, but for the singular values (i.e., the square roots of the eigenvalues 
of Q (fTT]) l. this requirement could be dropped. 



1.2.2. The Main Results of the Paper 

The first main result of our calculation, which we accomplish in section [2] 
by means of non-Hermitian planar diagrammatics, is the following simple 
equation for the non~holomorphic M-transform of P, 



n 



As already announced, the full information about the mean distribution of 
the eigenvalues of P can be retrieved in a straightforward way from this 
M-transform fOij) . (|A?2T]) . We see that ([9]) is a polynomial equation of 
order L. Moreover, since the RHS is a function of \z\^ only, so will the 
non-holomorphic M-transform be, 

Mp(z,z) = 9Jtp(|z|2). (10) 

Consequently, the spectral density will display rotational symmetry, in con- 
cord with the square case ([6]). Finally, remark that the linear scale of the 
distribution is still given by o" ([5]), i.e., z appears only as the combination 
zjo. 

The second main contribution of the article, worked out in section [3] 
with aid of the FRV multiplication law, is an equation obeyed by the (usual, 
Hermitian) M-transform of Q ([7]), yielding thus ()A.7p . (|A.6p the singular 
values squared of the product P, 



i (11) 



This is a polynomial equation of order (L-l-1). This equation for the product 
of square matrices has been derived in [12^ I13j. while for rectangular ones, 
without our knowledge, in [35] in the context of wireless telecommunication; 
we thus present our derivation of the latter result as well. 

As straightforward consequences of the main equations ^ and pT|) . 
we unravel the nature of the singularities of the mean spectral densities of 
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P and Q at zero. We show that their behaviors at the origin are solely 
determined by the number s of the rectangularity ratios ([2|) equal to 1, i.e., 

s = #{l = 1,2,..., L : Ni = Nl+i} = 1,2,..., L, (12) 

and are given by 

/op(A, A) ~ |Ar2(i-^), as A^O, (13) 

and 

Pq(A)~A"^, as A^O. (14) 

Remark that when all the constituent matrices are square, s = L, then (jl3p 
precisely reproduces the density in the entire domain of the spectrum 
not only close to the origin. 

Let us close with the following observation: (|lip happens to be remark- 
ably similar to @: Setting i?i = 1, which is a must if we wish to simul- 
taneously compute both the eigenvalues ([9]) and singular values (fTTjl . one 
notices that it is sufficient to replace on the LHS M-p{z,z) by Mq{z), and 
multiply by (Mq(z) + 1)/Mq(z), while on the RHS to replace |zp hy z — 
in order to proceed from ([9]) to (fTTj) . 

This unexpected striking relationship has led us to propose the follow- 
ing conjecture: If X is any non-Hermitian random matrix model whose 
mean spectrum possesses rotational symmetry, i.e., equivalently, whose non- 
holomorphic M-transform depends only on |zp (fTO]l . Mji_{z,z) = 9Jtx(|^P) 
— then, introducing the "rotationally-symmetric non-holomorphic A^-transform" 
as the functional inverse 

m^{m^{z)) = z, (15) 

its relation to the (usual) A^-transform of the Hermitian matrix X^X reads 

A^xtx(^) = ^^x(^). (16) 

One may intuitively expect the existence of such a link, because a rotationally- 
symmetric spectral distribution is effectively one-dimensional, depending 
only on the absolute value squared of its complex argument — and the 
eigenvalues of X"''X are precisely the absolute values squared of the eigen- 
values of X. From the practical point of view, typically, one of the two sides 
of ([16]) will be much easier to appropriate analytically, thus providing the 
other one for free. We postpone verification of this hypothesis for future 
work. 
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1.2.3. The Plan of the Paper 

The material in the article is divided into three levels of depth: 

• Experts in the field may remain just with the above ^1.2.21 which 
summarizes all of the results we have obtained in this work. 

• Intermediate-level readers may want moreover to consult sections [2] 
and [21 in which we compute the mean spectral densities of the product 
P ([H) of rectangular Gaussian matrices ([8]), as well as of Q = P^P ([7]), 
respectively. 

• Students are given furthermore appendices 1X1 and [Bl which introduce 
in a detailed and self-contained way non-Hermitian planar diagram- 
matics and free random variables calculus, respectively. 

Section HI concludes the paper and discusses some possible applications of 
its discoveries. 

2. The Eigenvalues of a Product of Rectangular Gaussian 

Random Matrices 

2.1. Derivation 

In this subsection, we present a derivation of the first main finding 
of this article, (HI), i.e., an L-th-order polynomial equation for the non- 
holomorphic M-transform — an object containing the full information 
about the mean spectrum of the (complex) eigenvalues (|A.24|) . (|A.2ip — 
of the (non-Hermitian) product P dH) of rectangular ([2]) — with the neces- 
sary constraint N^j^i = Ni — Gaussian random matrices dH). To this end, 
we employ an efficient technique of summing planar diagrams, called the 
"Dyson-Schwinger's equations," extended to the non-Hermitian sector, and 
specifically suited to working with products of random matrices (a reader 
not familiar with the subject is strongly advised to consult appendix \K\ for 
a didactic overview). 

2.1.1. The Linearization of the Product 

If one aims at computing the matrix-valued Green's function (|A.29p for 
a product P (HI) of random ensembles by means of planar diagrammatics, 
one notices that the problem becomes non-linear w.r.t. the constituent 
matrices. Opportunely, it is possible to linearize it by making use of the 
following trick: One trades the A'^i x A'^i matrix P for the A^tot. x Atot. — 
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where A'tot. = iVi + A''2 + . . . + A''l — matrix 



/ Ai 
A2 
























(17) 



with an L X L block structure, where the O's represent matrices of vari- 
ous appropriate dimensions entirely filled with zeros (we will also use the 
notation of Oat for the N x N zero matrix). 

Now, the non-holomorphic M-transforms of P and P are related in a 
simple way, 

Mp(t/;,uJ) = ^Mp(u;^,w;^). (18) 

We present the proof in ^A.2.31 In other words, solving the spectral prob- 
lem for P (jl7p is equivalent to solving the original model P ([T|). The 
former is already linear w.r.t. the constituent matrices, and permits the 
Dyson-Schwinger's approach to calculating Mp{w,W), from which the de- 
sired Mp{z,z) is obtained with help of (fTHj) . albeit for a price of enlarging 
the size of the random matrix from A'^i to iVtot.- 



2.1.2. The Dyson— Schwinger's Equations 

As demonstrated in ^A.2.2l the first step in writing the Dyson-Schwinger's 
equations^ is to know the propagators of the random matrix in question, 
namely, P, or more precisely, its "duplicated" version (|A.25p . 

P°-(oP°.)- <^'' 

It is a 2A'^tot. X 2A^tot. matrix, but we will think of it as having four blocks 
(as distinguished in (jl9p ). each being an L x L block matrix (see (jlTp ). 
We will denote the L x L block indices in these four blocks by Im (upper 
left corner), Im (upper right), Im (lower left), Im (lower right), each one 
covering the range 1, 2, . . . , L (compare (|A.36p ): for example, [P^J^r = A|. 
All the other featured matrices will inherit this structure, for instance. 
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[GD] 


12 


GD" 


IL 






11 


[GD] 


12 


[GD] 




21 




22 


gd; 


2L 






2T 




22 





[G°] 




G°] 




[G^ 


Tl 


[G^ 


T2 




21 




22 



[gd 



G^ 

gd 



G^ 

G^^ 



IL 
2L 



G^ 
GD 



IL 
2L 



11 
2T 



12 
22 



IL 
2L 



V [G°]zi [G^ 



\L2 



[G°]r. [G°]zT [G°]l^ ••• [G°]zz/ 

(20) 

and similarly for and 5]^. (We use w instead of z to comply with the 
notation in (jlSp . as well as disregard for simplicity all the subscripts and 
the symbols of dependence on w^w.) 

JtVe are interested in computing the non-holomorphic Green's function 
of P (TOO]) , i.e., 



Gp{w,w) 



TrG- = -i-5:TV[G°], = -3-j;iV,g,, (21) 



A^tot. ^tot- 



where it is useful to define the normalized traces 



(22) 



Hence, we should find the Qus. 

The 2-point correlation functions of the ensembles A/ are readily visible 
from their probability measures ([8]), 



At 



cd 



err 



--^Im^adhc 



(23) 



with all the rest equal to zero. In terms of P^, (|23p means that the only 
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non-zero propagators are 



pD" 


12 


pD" 




pD" 


23 


pD" 


32) 



pD" 




pD" 




LI 





IL 



lAfr 



LiVi, 



(24) 



where the tensor notation is a shortcut for what could be alternatively 
achieved using a double-index notation, 



pD" 




pD" 




(i,a),(/+l,A) 





err 



(i+i,B),{lb)/ y/NiNi+i 



5ab^AB, 



(25) 



where I = 1,2,...,L (and we use a cyclic identification convention here, 
L + 1 = 1), a,6 = 1,2,... ,A^i, and A,B = 1,2, . . . , Ni+i. 

Thus, we are now in position to write down the two Dyson-Schwinger's 
equations for P^. The first one, being in fact the definition of the self- 
energy matrix, is independent of the propagators and identical to (|A.38p 
and (|XT2]1 . 

(26) 



(WD-S°)"'. 



The second one is pictorially presented above formula ()A.13p , and the struc- 
ture of the propagators (|24p implies that the only non-zero blocks of the 
self-energy matrix read 



or 



i+i 



Ni 



-^i+iJ+i'^Ni, (27) 



=ai 



[^°] 



l-i 



^ i^'']— 1,1-1 = ^'lixl 1,1-1 l^n (28) 



for all / = 1, 2, . . . , L, with the cyclic convention = L, where the normalized 
traces ([22]) have been used. 

Results dsn, ([28]) mean that the four blocks of the matrix (W° - S°) 
are diagonal. 



W 



D 
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wIni 








-qiIatj 











w1n2 








—a2lN2 



















• -aL'i-NL 


-/3ll7Vi 








wIni 











—(32'i-N2 








w1n2 











• -Pl^Nl 








wInl 



(29) 



Such a matrix can be straightforwardly inverted: its four blocks remain 
diagonal, and read 











ai7ilAfi 











1(772 IaTj • 








0272 lA^a • 



















• aL7LlAfL 





























1^72 lAfa 






















(30) 



where for short, for all / = 1, 2, . . . , L, 



^ _ I |2 o I |2 / n2 V^^'-l^^'+V ^ 

- = M -aiPi = \w\ -{ai^iai) — Qi+^j^Qjzi^i_r (31) 

Substituting ([30]) . which is an implication of the second Dyson-Schwinger's 
equation, to the first one ()26p . we discover that the only non-zero blocks of 
the duplicated Green's function ([20]) are, for alH = 1, 2, . . . , L, 



[G%^ = /3ailN,, [G% = w^ilNr (32) 

Taking the normalized traces of both sides of every equality in ()32p . this 
leads to the final set of equations, 

Gil = wji, Gfi = am, Qji = /3ai, = w-fi- (33) 
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To summarize, the structure of equations ()33p is the following: The 
fourth one, as a general caveat, is the conjugate of the first one, so it is 
redundant. The second and third ones read 



2 



(34) 



S-ii = ^/-i^ ^^^irT,z-i7«- (35) 

We see that ([5T|) . ([M]) and form a closed set of 3L equations for 3L 
unknowns, G^j, Qj^ and 7; — these we will now ( ^2.1.3|) attempt to unfold. 
Once solved, in particular, when the 7/'s are found, the first expression in 
(f33|) yields Qii, i.e., as a consequence (f2T]) . the non-holomorphic Green^ 
function of P, and subsequently, the non-holomorphic M-transforms of P 
()A.24p . as well as of P (in the argument w^) (fT8|) . 

_ _ 1 ^ _ 1 ^ 

G^{w,w)=w— — ViVa/, «-e., M^{w,w) = — — Va^/W, 



i.e., Mp (i/;^ , lU^) = 1 ^ i?; w , (36) 

1=1 

where we have traded the 7;'s for a more convenient set of variables, 

m = \w\'^ji - 1, (37) 
and recall that the rect angularity ratios Ri = Ni/Ni 

2.1.3. Solving the Dyson— Schwinger's Equations 

Let us start from an observation that if we knew the 7i's, then equations 
([3^ and ([35]) would comprise a set of decoupled recurrence relations for Qj^ 
and Qji^ respectively. Assume this is the case, and iterate these recurrences 
down to / = 1, 

'\aia2...(Ti-iY ^/Rilil2...1i-i 

2 1 

Q-u = Q II {(yi(^2---m~i) ^=72... 7z- (39) 
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Actually, ([38]) and ([39]) hold true for all I = 1, 2, . . . , L (and not only 
/ > 2), where the / = 1 case is obtained by applying the cyclic convention 
= L. This cyclic boundary constraint for (|38p is 

(cricJ2 . . . aif (7172 • • • 7l) Gil = ^iT (40) 

(and analogously for ([39]) . with replaced by Gji)- One possibility now is 
that Qii = 0. This implies in turn ([38p that for all /, Qfi = 0, i.e., further- 
more dn]), 7« = l/ltL-P, or dSIl), /i; = 0, and therefore ([MD, Afp(z,z) = 0. 
This is thus the trivial holomorphic solution, valid outside of the mean 
eigenvalues' domain. We are however interested in the eigenvalues, i.e., in 
the non-holomorphic solution, so let us take Qij 7^ henceforth. The single 
equation ([^0]) reduces then to 

1 u^l' 

7i72---7L = ^, i-e-, (^1 + 1) (^2 + 1) • • • (/^L + 1) = — (41) 

with a defined in ([5]). 

Substitute now ([38]) and ([39]) into ([3T]). and go to the variables fii — 
after some simplifications, (|3ip acquires the form, for all / = 1, 2, . . . , L, 



_ H'^GijGii , . 

- w + l ' ^^^^ 



Take it first for / = 1, 



where we have used Ri = 1, and replace the RHS of ([I3]) appearing in ([^2]) 
by its LHS — obtaining thereby all the /i^'s, for / = 2, 3, . . . , L, as simply 
related to /ii, 

W = f • (44) 

Hence, if we manage to compute /xi, then all the /i;'s will be known as well; 
this is done by plugging ([44p into ([4ip . which then becomes a polynomial 
equation of order L for /xi. 



This completes the solution of the fundamental set of equations ([HI]) . ([MD, 
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As mentioned, the knowledge of all the fiis yields the desired non- 

holomorphic M-transform of P; indeed, substituting into (p6|) . we 
simply get 

Mf {w^ ,w^) = ^ii. (46) 

In other words, it remains to change the complex argument from w to 
z = in order to see that Mp(z, z) obeys the L-th-order polynomial equa- 
tion 

which is precisely the first main result of our article, Q. 

The last point is to determine the domain of validity of the non-holomorphic 
solution (j47|) . i.e., the domain where on average the eigenvalues of P lie. We 
know ()A.23p that on the boundary of this domain, the non-holomorphic and 
holomorphic solutions must agree; plugging thus the latter (Mp(z,z) = 0) 
into (HZI, we obtain the equation of this borderline, 

\z\=a. (48) 

It means that the eigenvalues of P are scattered on average within a centered 
circle of radius a, with the density stemming from (j47p by virtue of ()A.24p . 
(TOTD . 



2.2. Comments 

2.2.1. Example: L = 2 

Since our main equation ()47p is polynomial of order L, in the special 
case of L = 2 it will be easily solvable. Indeed, the non-holomorphic M- 
transform reads 



Mp(z,^) = ^-l-i?+y(l-i?)2 + 4i?^j , (49) 

where we call R = R2 = N2/N1 (the only non-trivial rectangularity ratio 
here), and where we have picked one solution out of the two roots of the 
quadratic equation ([T7|) in such a way that it satisfies M-p{z,'z)\\^\^^ = 0, 
i.e., the matching condition with the holomorphic solution on the borderline 
(|48p . As a result, we immediately arrive at the non-holomorphic Green's 
function (lAlMl) . 



Gp(z,z) = ^ fl-i?+y(l-i?)2 + 4i?^ ) . (50) 
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Fig. 1. Numerical verification of tfie theoretical formula ((52|) for (the radial part 
((57|) of) the mean spectral density /Op(z,z) of the product P of L = 2 rectangular 
Gaussian random matrices, as well as the finite-size correction (j58p . A numerical 
histogram (the black line) versus the theoretical prediction (|52p . supplemented 
with the finite-size smoothing (155)) (the red plot), for A^i — 100 and N2 — 200 
(i.e., R — R2 ~ 2), and for 10^ Monte-Carlo iterations {i.e., the histogram is made 
of 10^) eigenvalues. The adjustable parameter q ([55| is fitted to be q « 1.14. 




Fig. 2. An analogous graph to figure[l] this time with Ni = 100 and N2 = 150 {i.e., 
R = 1.5). We find q w 1.08 here. 
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Fig. 3. An analysis of the finite-size effects: Numerical histograms for Ni — 50, 
N2 = 100 (black), iVi = 100, N2 = 200 (dashed red), iVi = 200, N2 = 400 (dotted 
blue), i.e., with the same rectangularity ratio R — 2, but increasing matrix di- 
mensions. We observe how these plots approach the green line of the theoretical 
formula (|52p for the density in the thermodynamic limit. 
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Fig. 4. Numerical histograms for the matrix sizes of Ni = 100, A'2 = 200 {i.e., 
R = 2; black) and iVi = 200, N2 = 100 {i.e., R = 1/2; red). Due to the presence 
of the zero modes (not displayed in the picture), the latter is half of the former. 



The average density of the eigenvalues stems from (|5Up by taking the 
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Fig. 5. Analogous graphs to figure [31 but for L — 3 (up) and L — 4 (down). 



derivative w.r.t. z (|A.2ip . however, one has to be cautious in the vicinity 
of zero in order to properly take into account the zero modes. Let us first 
expand ()50p near z = for the purpose of making its behavior clearly visible. 



Gp(z,z) ~ ^ +regular terms, where /=| q ^' for > l' 

" (51) 

Applying the derivative (l/vr)^^ to this singular term produces the complex 
Dirac's delta at the origin, fS^^\z,z) (see the discussion after ()A.53p ) . Away 
from z = 0, however, this term is irrelevant, and only the square root in 
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([50|) contributes to the density — i.e., altogether, 

for l^l < fj, 

(52) 

for \z\ > a. 

As for any L, this function (|52p possesses rotational symmetry. One 
may check that it is correctly normalized, Jq dr2'Krp-p{z,'z)\\^\^.^ = 1, the 
delta function playing a vital role in this agreement. Setting R = 1 gives 
/9p(z,z) = l/(27r(T|2;|), for \z\ < a, in concord with the finding ([6]) of |14j . 



R 



PP{Z,Z) 



(l-_R,)2+4iJ^ 



+ f6('Hz,z), 



0, 



2.2.2. The Singular Behavior of the Mean Spectral Density at Zero 

Let us now examine one property of the mean spectral density of P, 
namely, its behavior at the origin of the complex plane. We will show that 
the density is singular there, and that the rate of this explosion is governed 
by the number s = 1,2, ... ,L ()12p of the rectangularity ratios Ri equal to 
1. 

Indeed, we begin from recasting the main formula ()47p in the language of 
the non-holomorphic Green's function ()A.24p . and explicitly distinguishing 
the terms with Ri = 1 from those with Ri ^ 1, 

ie{l....,L}: ^ ' 



Consider now how ()53p behaves in the limit of z — )• 0. Let us suppose 
that as z — >• 0, so also zG-p{z,z) — ?■ 0; we will justify this assumption a 
posteriori. Consequently, the entire product over I such that i?; 7^ 1 in ([53]) 
tends to a non-zero constant, and therefore, the singular behavior of the 
non-holomorphic Green's function reads 

z^Gp{z,zy '-^ zz, i.e., Gp{z,z)'~~^z^'^z^, as z — )■ 0. (54) 

(We confirm that our previous assumption holds, zG-p{z,z) ~ |zp/*.) Tak- 
ing the derivative w.r.t. z of both sides of ()54p . the singular behavior of the 
density (|A.2ip is finally found, 

pp{z,z) \z\~'^(^~^\ as z ^ 0, (55) 

as anticipated in (|13p . (Actually, for s = 1, and only for this value, there is 
no singularity, pp{z,'z) — >• const.) 
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2.2.3. The Behavior of the Mean Spectral Density at the Borderline 

Let us also focus for a while on the behavior of the density pp(z,z) at 
the opposite end of the spectrum, namely, at the centered circle of radius a, 
constituting the borderline of the average eigenvalues' domain. It has there 
a finite value of 



1 11 

pp{z,z) = — 2^h, where ^ = (56) 

(Note that it complies with ([52]) and ([6]).) To see this, recall that for the 
rotationally-symmetric non-holomorphic M-transform, Mp{z,z) = Wlp{\z\'^) 
the density is related to it as /9p(z, z) = (l/7r)a|^|29Jtp(|z|2). Taking now the 
derivative w.r.t. Mp{z,z) of both sides of our main equation ()47p . and set- 
ting Mp{z,z) = {i.e., the value at the boundary), we readily arrive at 

Therefore in the considered thermodynamic limit ([5]) the density has 
a jump from the value (|56|) to zero as one crosses the frontier \z\ = a. 
However, for finite sizes of the random matrices in question, this step gets 
smoothed out. Imitating the idea of [SSl [371 [38] , we propose the following 
model for this finite-A^ (where by N we denote the order of magnitude of 
the dimensions of the matrices, say N = Ni) smooth crossover from within 
the eigenvalues' circle to its outside: Since the mean spectral density of P 
is rotationally symmetric, it is sufficient to consider its radial part, 

pj?'i-(r)^2^rpp(z,z)||,|=, . (57) 

Subsequently, we introduce the following "effective" radial density, which is 
supposed to properly incorporate the finite-A^ behavior at the borderline, 

pf-{T) ^ /5jp^^ (r)ierfc (g(r - l)^/iv) , (58) 

where the complementary error function is erfc(x) = {2/ y/ir) dt exp(— t^), 
while g is a free parameter, whose value is to be adjusted by fitting. We 
numerically verify this hypothesis (j58p in the next paragraph. 



2.2.4. Numerical Confirmation 

We have tested the main equation (j47p . as well as the finite-size ansatz 
(|58p , quite extensively, obtaining excellent confirmation in all the considered 
situations, see figures [U [21 [3 [D and [SI 
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3. The Singular Values of a Product of Rectangular Gaussian 

Random Matrices 

3.1. Derivation 

This subsection is devoted to deriving tiie second main formula of this 
article, (jlip . i.e., an (L+l)-th-order polynomial equation obeyed by the M- 
transform, a quantity equivalent to the average density of the eigenvalues, of 
the Hermitian matrix Q = P^^P ([7]), with P the product ([T]) of rectangular 
^ Gaussian random matrices ([5]) . The underlying idea will be to rewrite Q 
as a product of some Hermitian matrices, to which the FRV (appendix iB]) 
multiplication law ()B.20p will be harnessed. 

Let us commence from defining, for any / = 1, 2, . . . , L, a square NiJ^i x NiJ^i 
matrix 

= (Ai A2 . . . A,_i AO^ (Ai A2 . . . Az_i A/) = 

= a\a\_^ . . . A^aJ A1A2 . . . A,_iA;, (59) 

being a generalization of Q which includes only the first I random matrices, 
as well as a square Ni x Ni matrix, which difi^ers from Q; only by the position 
of the last matrix in the string, i.e., A;, which is now placed as the first 
matrix in the string, 

Q; = AiA\a\_^ . . . At a1 A1A2 . . . A,_i = (a^aI) Q,_i. (60) 

We are interested in the eigenvalues of the Hermitian matrix Q = Ql- 

The orders of the terms in the two above products (f59|) . ([60]) are related 
to each other^by a cyclic shift, therefore, for any integer n > 1, there will be 
TrQ" = TrQ". Hence, the M-transforms (jA.Sp of the two above random 
matrices are related by 

n>l 7i>l 

(61) 

Inverting ()6ip functionally, we obtain a relationship between the respective 
A^-transforms (jB.lOp . 

The reason for introducing the auxiliary matrix (|60p is that it is a 
product of two free matrices, A;A| and Qi_i. Using the FRV multiplication 
law (|B.20p . we can thus write, for all / = 2, 3, . . . , L, 

^Q.(^) = JTT^A'^I^"^^^'-^"^- ^^^^ 
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From these two equations (|62p , ()63p , we now eliminate the A^-transform 
of the auxihary Q/, leaving us with the following recurrence relation for the 
A^-transform of Q;, 



^^Q.(-) = ^^^^A.At(^-)^Q.-.(^-), for / = 2,3,...,L, 

(64) 



with the initial condition, 



A'Q,M = A'«.(f^)='V...,(|^). (85) 

which stems from (j62|) and the fact that Qi = AiA|. The solution of this 
recurrence is then readily found to be 



(z + R2) iz + R3)...{z + Rl) 



■'V....(j^)'V.,..(i)...A'.,4(^). W 

It remains now to find the A^-transforms of the random matrices A;A|. 
They are examples of the so-called "Wishart ensembles," and the problem 
of computing their A^-transforms, with the precise normalization of the 
probability measures of the A^'s which we are employing ([8]), has first been 
solved in [39]: Expressions (1.8), (2.8), (2.13), (2.14) of this article yield 
the Green's function of A^aJ, which immediately leads to the pertinent 
A^-transform, 



(-+i)U/ie7-+ 



^A,At(^) = -I ^ ^ -■ (67) 



Substituting (|67|) into (|66|) . one finally arrives at the desired formula for the 
A^'-transform of Q = Ql, 

^v«(.) = + 1) (^ + 1) (^ + 1) . . . (^ + 1) . m 

with fj defined in ([5]). In other words, the corresponding M-transform 
M(^{z) satisfies the following polynomial equation of order (L + 1), 

(69) 
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or in the case of N^j^i = Ni {i.e., Ri = 1, required when one wishes for P 
to have eigenvalues too), 



M, 



This completes our derivation of ([TT 



(70) 
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Fig. 6. Numerical verification of the theoretical formula (j69p for the mean spectral 
density Pq(A) of the random matrix Q = P^P ([7]). Everywhere Nl+i — Ni — 50. 
The number of Monte-Carlo iterations is 20,000, i.e., all the histograms are gen- 
erated from fO^ eigenvalues. Here L — 2, and the matrix sizes are chosen to be 
Ni = 50, N2 = 150. 



3.2. Comments 

3.2.1. The Singular Behavior of the Mean Spectral Density at Zero 

Following an analogous line of reasoning as in ^2.2.21 we may unravel 
the singular behavior at zero of the mean spectral density of Q stemming 
from the main formula ()69p . Again, we rewrite this equation in terms of the 
(usual, holomorphic) Green's function ()A.7p . and separate the terms with 
Ri = 1 (the number of such rectangularity ratios is called s ()12p ) from those 
with Rij^l, 

zGci(z) _ \\ z 
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Fig. 8. The same as figure [SI but with L = 4, and the matrix sizes are A^i = 50, 
TVs = 100, N-i = 150, 7V4 = 200. 



In the limit z — t- 0, the relevant terms in (|7ip will thus read 

z^+1Gq(z)'+i ~ i.e., Gq(z)~z"^, as z ^ 0, (72) 
which immediately leads to the aforementioned result (|14|) . 
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3.2.2. Numerical Confirmation 

We have performed several numerical tests of the main formula ()69p . in 
all cases obtaining astonishing agreement, see figures El [7] and [51 

4. Conclusions 

4-1. Summary 

The main contribution of this article is given by the two equations ([9|) 
and (|lip for the M-transforms, i.e., objects conveniently encoding the aver- 
age spectral densities, of the product P = A1A2 . . . A.i ^ of an arbitrary 
number L of independent rectangular (|2]) Gaussian random matrices ([8]), 
as well as of the matrix Q = P^^P (jT]), whose eigenvalues are the squared 
singular values of P, respectively. Both these equations are polynomial, and 
have orders L and (L + 1), respectively, so in general they may be solved 
only numerically; however, some properties of the mean spectral densities 
can still be retrieved analytically, such as their singular behavior at zero 

Furthermore, (jH) and (jlip are very similar to each other — the fact which 
directed us to put forward a conjecture that the same resemblance is a 
feature shared by all non-Hermitian ensembles X possessing rotationally- 
symmetric average distribution of the eigenvalues. For such models, the 
non-holomorphic M-transform Mx.{z,'z) ()A.24p is a function of the real 
argument |zp pop . thereby allowing functional inversion, and hence, a defi- 
nition of the "rotationally-symmetric non-holomorphic iV-transform" (|15p 
— even though for general non-Hermitian random matrices a construction 
of a "non-holomorphic A^-transform" remains thus far unknown. This new 
A^-transform is then conjectured to be in a simple relationship ()16p to the 
(usual) A^-transform of the Hermitian ensemble X^X. In a typical situation, 
the latter will be solvable much more easily than the former, owing to the 
plethora of tools devised in the Hermitian world, albeit the opposite may 
be true as well. This is indeed the case here — our derivation of Q, based 
on non-Hermitian planar diagrammatics and Dyson-Schwinger's equations, 
is much more involved than a simple application of the FRV multiplication 
rule leading to ([TT]) — and consequently, the hypothesis would provide 
a means of avoiding the complicated diagrammatics. To the best of our 
knowledge, this would be the first use of the free random variables calcu- 
lus to computing the mean spectral density of a non-Hermitian product of 
random matrices. 

We have also suggested a model for a finite-size behavior of the density 
of P near the borderline of the eigenvalues' domain (|58p . taking after the 
work of |36t [571 138| . It performs outstandingly well when checked against 
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numerical simulations. 



Possible Applications 

Let us now succinctly sketch some possible applications of these results 
to wireless telecommunication, quantum entanglement and finance. 

4.2.1. Wireless Telecommunication 

Information theory for wireless telecommunication has been intensively 
developed in the past decade, after it had been realized that the infor- 
mation rate can be increased by an introduction of multiple antenna chan- 
nels, known as the "multiple-input, multiple-output" (MIMO) transmission 
links. The MIMO capacity for Gaussian channels has been calculated in the 
pioneering work |40j . triggering large activity in the field. Immediately, it 
became clear that an appropriate language and methods to address this 
type of problems are provided by random matrix theory (consult [6] for a 
review) . 

The MIMO capacity reads [lO], 



Capacity C = (log^ Det (^Ijv,,,. + ^ ACAt^ ^ , (73) 

where A'rec. is the number of receivers, Ntr. of transmitters, and SNR the 
signal-to-noise ratio. The output signals y are calculated from the input x 

as 

/SNR^^ ^ 
^"V^ (74) 

where A is the response matrix for a given frequency, r/ is a standardized 
multivariate white noise, while C is a covariance matrix for the input signals, 

i.e., [C]ab = {XaXb). 

In the simplest case, one assumes that the input signals are uncorre- 
lated, i.e., C = iTVtr > ^iid that A is a random matrix built of IID cen- 
tered Gaussian entries. This corresponds to a random situation, when 
one has no information about the signal propagation. Then, the asymp- 
totic (A'tr. — >• oo) mutual information per channel can be derived to be 
H = j dA/?Q(A) log(l + SNRA), where /Oq(A) is the limiting mean eigenvalue 
density of the Wishart matrix Q = A"^A [6j. 

The model considered in our paper can be applied to a more complex 
case of signals traveling over L consecutive MIMO links [35]: It is first 
sent from Ni transmitters via a MIMO link to receivers, which then 
re-transmit it via a new MIMO link to the subsequent A^3 receivers, etc. 
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Clearly, the capacity will depend on these numbers of intermediate re- 
transmitters; in particular, if any of the A'^^'s is small, the capacity will 
be reduced. The output now reads 



and therefore, the effective propagation is given by the matrix P = Aj^ . . . A2A1, 
whereas the mutual information per channel hy m = J dA/9Q(A) log(l + SNRA), 
with Q = P^^P. Now, it is precisely our second main result (jlip which can 
be exploited to find the relevant density Pq(A). 

Let us finally just mention that one could imagine a more general sit- 
uation, where MIMO links form a directed network — each directed link 
Im representing a single MIMO channel between A'^; transmitters and Nm 
receivers. (The discussed case ([75]) corresponds to a linear graph, 1 — )• 2 — )• 



4.2.2. Quantum Entanglement 

A complex directed network of MIMO links is somewhat similar to 
the structures appearing in the context of quantum entanglement. There, 
one considers graphs whose edges describe bi-partite maximally entangled 
states, while vertices — the couplings between subsystems residing at the 
same vertex • In the simplest case of a graph consisting of a single link, 
it is just a bi-partite entangled state. The corresponding density matrix 
for a bi-partite subsystem is given by Q = A^^A, where A is a rectangu- 
lar Ni X N2 matrix defining a pure state, IV') = X^^i Z^b^i[A]afe|aa) <^ \Pb), 
being a combination of the basis states in the subsystem, \aa) and |/3;,) (see 
for instance [42]). One can easily find that linear graphs with additional 
loops at the end vertices correspond. The density matrix for the subsystem 
sitting in the end vertex is given by Q = P^P, where P = A1A2 . . . |41| . 
If all the subsystems are of the same size, the average spectral distributions 
of Q are known \T2[ [T3\ as the "Fuss-Catalan family" [43j; they can be 
obtained from ()lip by setting all the i?;'s to 1. However, if the subsystems 
have different sizes, one needs to apply our general formula (jlip . 

4.2.3. Finance 

Another area of applications is financial engineering. Products of rect- 
angular random matrices may be used in calculations of lagged correlation 
functions, which play very important roles in risk management and portfolio 
theory. This issue will be discussed in a forthcoming publication. 




(75) 
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Appendix A 

The Diagrammatic Approach to Solving Matrix Models 



Appendix A.l Hermitian Random Matrices 
Appendix A. 1.1 The Green's Function and A<f— Transform 

The most basic question one asks in RMT about an A'^ x A'^ — where we 
take the thermodynamic limit of A'^ — > oo — Hermitian random matrix H, 
endowed with some probability measure d//(H), is to compute the average 
distribution of its (real) eigenvalues Ai, A2, . . . , Aat, defined as 

1 ^ 

Ph(A)^-5^(<5(A-A„)), (A.1) 

a=l 

where the average map (...) = J'(. . .)d/x(H), while (5(A) denotes the real 
Dirac's delta function. 

However, it proves more convenient to work with another equivalent 
object — the "Green's function" (a.k.a. "resolvent"). First, one introduces 
the N X N "Green's function matrix," 

Gh(z) = <(Z-H)-^), where Z = z1n, (A.2) 

where z is a complex argument, and In represents the N x N unit matrix. 
Then, the Green's function is defined as its normalized trace, 

1 1 ^ / 1 \ 

G„(.).-T.G„(.) = -!:( — ). (A.3) 

a=l ^ 

As clearly seen in this definition, Gh(-z) is, for any finite N, a meromorphic 
function, with poles located at the average spectrum. When A^ tends to 
infinity, these poles coalesce into continuous intervals on the real axis, and 
Gh(-z) turns into a function holomorphic everywhere on the complex plane 
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excluding the cuts formed by these intervals. Knowing the mean density of 
the eigenvalues in these cuts allows one to reproduce the Green's function, 

Gyi{z) = [ dApH(A)^ (A.4) 

(in other words, the Green's function is the Stjeltes' transform of the den- 
sity), and conversely, a standard representation of the real Dirac's delta, 

d(X) = -- lim Im-^ (A.5) 

^ ^ 7re^o+ A + ie ^ ^ 

implies that the mean spectral density is recovered when one approaches 
with the argument of the Green's function the eigenvalues' cuts along the 
imaginary direction, 

/oh(A) = -- lim ImGnfA + ie). (A.6) 

We thus see that (jA.ip and ()A.3p carry identical information. 

Instead of the Green's function, one often prefers yet another quantity, 
the "M-transform," simply related to the former, 

Mh{z) = zGn{z) - 1. (A.7) 

Its meaning is deduced from considering its large-z expansion — being 
holomorphic everywhere except the cuts on the real line, it permits a power- 
series expansion around z = oo, 

Mn{z) = where Mh,„ = ^ Tr (H") = / dApH(A)A". 

n>l ^ J cuts 

(A.8) 

The coefficients are precisely the moments of the probability distribution 
in question, and so, Mii(z) is also known as the "moments' generating 
function." 



Appendix A. 1.2 The Green's Function for the GUE from Planar Diagrams 

There are multiple techniques developed for the purpose of calculating 
the Green's function of a Hermitian random matrix H, and we will now 
briefly describe the method of representing it as a sum of planar fat di- 
agrams, which is then performed using the so-called "Dyson-Schwinger's 
equations." We will not aim at a comprehensive presentation, but only 
show how the approach works on a simplest example of a Hermitian model 
— the "Gaussian Unitary Ensemble" (GUE), with standard deviation a, 

d/i (H) oc e"^ DH. (A.9) 
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(It means that the real elements on the diagonal are IID with the same 
variance cr^/iV, while the real and imaginary parts of the entries below the 
diagonal are IID with the variance a'^/2N.) 

One begins from expanding the Green's function matrix ()A.2p into the 
power-series around infinite z, 

Guiz) = + (Z^^HZ^^HZ^^) + (Z^^HZ^^HZ^^HZ^^HZ^^) + . . . , 

(A.IO) 

where we have preserved only even moments, as the probability measure 
()A.9p is even, as well as kept the order of the Z~^'s and H's intact, even 
though they commute. Since Z is independent of H, and can thus be taken 
outside of the average, one is left at every order with an expectation of the 
form ([H]aia2[H]a3a4 • • • [H]a2„-ia2„ ) • Owing to Wick's theorem, any such 
n-point correlation function can be expressed — by making all the possi- 
ble contractions — through the 2-point correlation function (propagator), 
which for the GUE measure ()A.9p reads 

2 

{[UUiUU) = ^Saddbc (A.ll) 

Consequently, the structure of this series (|A.10p lends itself to a graphical 
representation, in which every factor of [Z^^J^fe is depicted as a straight line 
connecting the matrix indices a and b, while every propagator ([H]Q5[H]c(i) 
as a double arc connecting two pairs of indices, ab and cd, 




The first three orders of (|A.10p are thus drawn as 




I I I I 

([H]ciC2 [H]c3C4 [H]c5C6 [HjcTCg) 
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(There is summation over all the internal indices.) At the first order, we just 
have one horizontal line corresponding to [Z~^]af,; at the second one, there 
are three Z~^'s (i.e., three horizontal lines) and one propagator {i.e., one 
double arc); while the third order contributes five Z~^'s (i.e., five horizon- 
tal lines), and a 4-point correlation function, which permits three possible 
contractions, giving rise to three "rainbow" diagrams, down to products of 
two propagators (i.e., two double arcs), etc. 

Clearly, complete evaluation of [GH(-z)]afe5 for any finite matrix size A^, 
would require summing up the numerical values stemming from all such 
connected rainbow diagrams with external indices a and h. This is where 
the thermodynamic limit of A^ — )• oo comes into play — as is well-known 
since the work of 't Hooft [H], in this limit, only planar graphs contribute 
to the Green's function, while all the non-planar ones are suppressed with 
the factor of l/N"^^, where h is the genus of the surface on which a given 
diagram can be drawn without crossing of lines. For instance, the last graph 
pictured above can only be drawn without a crossing on a torus (genus 1), 
so it may be safely neglected at large A^. 

Now, all the planar diagrams can be summed up by exploiting a cer- 
tain technique known from quantum field theory, which we will now sketch. 
First, one considers a subset of the rainbow diagrams, referred to as "one- 
line-irreducible" (ILI) ones, that cannot be split into two disconnected 
pieces by cutting a single horizontal line. Let J^uiz), called the "self-energy 
matrix" (it is an A x A matrix), denote their generating function. 
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One now perceives that any planar rainbow diagram can be constructed 
by putting together horizontal lines and ILI graphs, i.e., a relationship 
between the full and ILI generating functions pictorially reads 





+ : 





+ 



I.e., 



Gn{z) = Z-i+Z-15]h(^)Z-i+Z-1Sh(^)Z-1I]h(^)Z-i+. . . = (Z - ^H{z)y 

(A.12) 

This is the first Dyson-Schwinger's equation. Remark that it looks very 
similar to (|A.10p . (|A.2p . except for the now-absent averaging; thus, we 
may say that the self-energy matrix represents an effective matrix H which 
allows one to remove the averaging operation from (|A.10p . ()A.2p . 

Another observation is that if one takes an arbitrary planar rainbow 
diagram and adds an external double arc to it, one obtains a ILI graph; 
and conversely, any ILI graph has a form of a double arc encircling some 
diagram. Hence, the generating functions in question satisfy 




which is 



C1,C2 = 1 ^ ^ 



RectGGProd printed on January 15, 2013 



33 



i.e., T.^{z) = a''Gu{z)lN. (A.13) 

This is the second Dyson-Schwinger's equation. Since it reveals that the 
self-energy matrix is proportional to the unit matrix, one may take the 
normalized trace of (|A.13p and rewrite it as 

T.u{z) = a^Gu{z), (A.14) 

where, analogously to (|A.3p . the "self-energy," 

SH(^) = ^Tr5]H(^). (A.15) 

Finally, substituting Sh(2) from the second Dyson-Schwinger's equa- 
tion (jA.lSp into the first one (|A.12p . and applying the normalized trace to 
both sides — yields a quadratic equation for the Green's function of the 
GUE ensemble (|X9]l . 

Gii{z)[z-<y^Gii{z)) =1, i.e., Gh(^) = ^ - VF^2^\/FT2^) , 

(A.16) 

where we have selected the correct solution out of the two by referring to the 
condition Gh(-z) ~ for z — )• oo (equivalent to the normalization of the 
density, /^^^^ dApH(A) = 1, see (|A.4p ): also, the square roots are in the prin- 
cipal branches. The average eigenvalue density is found from this Green's 
function using ()A.6p . and is the famous "Wigner's semi-circle distribution," 

r\/4cT2 - A2, for \e[-2a,2a], 



^^^^ \ 0, for A gM\ [-2a,2a]. ^ ^ ^ 

On this pedagogical example, we have shown how planar diagrammatics 
and Dyson-Schwinger's equations can be used to solve Hermitian matrix 
models in the thermodynamic limit. Let us now see whether non-Hermitian 
ensembles could be appropriated along similar lines. 

Appendix A. 2 Non~Hermitian Random Matrices 
Appendix A. 2.1 The Non— Holomorphic Green's Function and A^— Transform 

The chief difference between Hermitian and non-Hermitian random ma- 
trices is that the eigenvalues Ai, A2, • • • , Aat of the latter (call it X) are 
complex — in the thermodynamic limit of — >• 00 merging into some 
two-dimensional "domains" T) on the complex plane — with their average 
distribution given by the complex Dirac's delta, 

;,x(A, A) = - ^ (5(2) (A - A„ A - a:)) . (A.18) 

a=l 
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Thus, the concepts developed in ^A.l.ll and relying on the reality of the 
eigenvalues, the Green's function and M-transform, lose their meaning. 

The bottom line is that the complex Dirac's delta has distinct properties 
from its real counterpart. For instance, its useful representation will now 
be 

5(2)(A,A) = ihm ^ = -4=lim ,^„^ (A.19) 

as contrasted with ()A.5p . Mimicking this form, one is led to define the 
"non-holomorphic Green's function" as 

1 ^ / z-Y 
G-j^iz,^) = lim lim " 



a=l * I ui I 

= lim lim — Tr / ^Ijv - Xt \ 

e^QN^ooN \(zl7V-X)(zljv-Xt) + e2ljv/ ^ ' 

(we have ambiguously written here the matrix inversion as a fraction — we 
assume henceforth A/B = AB~^, i.e., multiplication by the inverse matrix 
from the right), since then the density (jA.lSp is just proportional to its 
derivative w.r.t. z, 

P^{z,z) = -^G^{z,z), for zGP. (A.21) 
vr cTz 

Notice the order of the limits in ()A.20p — for finite A^, one could immediately 
set the regulator e to zero everywhere except a finite set of points Aa, and 
the non-holomorphic Green's function would reduce to the usual one ()A.3P : 
only for infinite N and inside the eigenvalues' domains T>, the limit e — t- 
yields a non-trivial object, carrying information about the mean eigenvalue 
density. Outside P, one is always left with the usual Green's function ()A.3p . 

G^{z,z) = G^{z), for z^V, (A.22) 

which however does not tell this time anything about the eigenvalues. Ac- 
tually, this formula holds on the boundary of V as well, 

Gx(z,z) = Gx(^), for zedV, (A.23) 

which we will see to be precisely the equation fulfilled by the coordinates of 
the boundary of V. 

Similarly as in the Hermitian case, one finds it useful to consider a 
related quantity, named the "non-holomorphic M-transform," and defined 
analogously to ()A.7p . 

M^{z,z) = zG^{z,z) -I. (A.24) 
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Again, outside of T) it simplifies to the usual M-transform Myi{z), which is 
however obviously incapable of capturing the distribution of the eigenvalues 
of X. 

Important as it is, an essential hindrance in an effective usage of the 
non-holomorphic Green's function ()A.20p is the quadratic (in X) structure 
of its denominator. As a response to that, the following linearization trick 
has been proposed: Introduce a 2A^ x 2A^ matrix 

Gg(2,J) E limjim ((Z? - X°)-') , (A.25) 



where 



iel^ zIn )' \Qn Xt 

The superscript "D" comes from "duplication." Considering the four N x N 
blocks of this matrix, distinguished by the superscripts zz, z'z, zz and Yz 
(and omitting the symbols of the functional dependence on z,z, which we 
will often do in similar expressions), 

^x / ' (A.27) 



as well as its normalized "block-trace," 



, A B \ / TrA TrB , 



C I) J ^ \ TrC TrD 

i.e., a 2 X 2 matrix, referred to as the "matrix-valued Green's function," 
consisting of the normalized traces of the four blocks, 

one recognizes that the upper left corner of the matrix-valued Green's func- 
tion is exactly equal to the non-holomorphic Green's function ()A.20p . 

g^{z,z) = G^{z,z). (A.30) 

In this way, the relevant non-holomorphic Green's function has been en- 
coded as a part of another object, having linear structure in X, and looking 
very similar to the usual Green's function (compare ()A.25p with ()A.2p ). the 
price to pay for that being "duplication" of the matrices. This will enable 
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us to use techniques such as planar diagrammatics to compute the matrix- 
valued Green's function of a non-Hermitian ensemble X in a manner parallel 
to how one does it in the Hermitian sector ( ^A.1.2p : see below ( ^A.2.2p . 

Let us close with the following remark: We have seen that the zz- 
component of the matrix-valued Green's function is equal to the non- 
holomorphic Green's function ()A.30p . and thus, inside T) reproduces the 
mean eigenvalue density (|A.2ip . while outside P reduces to the usual Green's 
function (jA.22p . The zz-component is just its complex conjugate, = G^, 
so it carries no additional information. One might wonder about the off- 
diagonal elements, 

(z,z) = lim lim — Tr / ttv^ ^ 5 V (A.31) 

g¥iz,z) = lim lim — Tr ( n ) • (A.32) 

^xv , ; e^ojv^ooiV \(zl7v-X)(zliv-Xt) + e2l7v/ ^ ' 

Since they are equal to each other and purely imaginary, \Q^{z,'z) = \Q^{z,'z\ 
it is convenient to consider their negated product, which is a non-negative 
real number, 

C^{z,z) = -gii{z,z)gii{z,z) e M+ U {0}. (A.33) 

Clearly, it vanishes outside of the domains of the eigenvalues V, as the 
regulator may be safely zeroed there, 

Cx(z,z) = 0, for z^V, (A.34) 

whereas inside D, it acquires non-trivial (strictly positive) values. Conse- 
quently, it contains all the information about the boundary of P — if one 
finds C:s.iz,z) inside V, and subsequently sets it to (|A.34p . one obtains an 
equation obeyed by the coordinates (x,y), where z = x + iy, of dV. (Note 
that this is equivalent to ()A.23p .) This being our primary usage of Cx.{z,z), 
let us nevertheless mention that it has yet another content — it describes 
certain statistical features of the left and right eigenvectors of X [45j . 



Appendix A. 2. 2 The Non— Holomorphic Green's Function for the Girko— Ginibre 
Ensemble from Planar Diagrams 

Since the matrix-valued and usual Green's functions have analogous 
forms, the difference being the duplicated structure in the latter case, planar 
diagrammatics and the Dyson-Schwinger's equations will work in a similar 
way. Instead of giving a general introduction, let us present the method on 
a simplest example — the Girko-Ginibre matrix model A (with standard 
deviation a), defined in 1^. 
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The propagators are immediately read off from tliis measure 



\ab 



At" 


cd) ( 


At" 










ab 



\cd 



([A]a6[A] 



cd) 



At 



ab 



At 



cd 



0. 



(A.35) 



However, in order to employ the diagrammatic tools, one should know the 
propagators of the duplicated matrix A^ ()A.25p . Denoting the indices in 
its four N X N blocks by ab, ab, ab, ab, respectively, mimicking the labeling 
of zz, z'z, zz, in (|A.27p . one sees that the only non-zero propagators read 

([A°]., [A°],,> = ^5-,5,„ <[A°],, [A^]J = ^6,,6,,. (A.36) 



N 



This is all there is required in order to write down the Dyson-Schwinger's 
equations. The self-energy matrix will be duplicated. 



•^zz 



(A.37) 



(for simplicity, we remove the subscript "A" in this paragraph). The first 
Dyson-Schwinger's equation will clearly have an identical form to ()A.12p . 



(A.38) 



where TP = Z^g ()A.26p — the regulator e may henceforth be set to zero; 
we will see that the Dyson-Schwinger's equations by themselves take care of 
properly regulating the result inside the mean eigenvalues' domains T). The 
second Dyson-Schwinger's equation will also pictorially look like before (see 
the diagrammatic equation above ()A.13P ). but the present structure of the 
propagators ()A.36p implies that the four blocks of that formula now read 



J ah 



N 

Cl,C2 = l 



[A°].. [A°] 



C2b 



a 



TtG' 



[^\b = 



2rzz : 



i.e., altogether. 



Otv 



TV 



Otv 



^ab = -''Q^'^Kb^ 

(A.39) 
(A.40) 

(A.41) 
(A.42) 
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Substituting the expression ()A.42p for the duphcated self-energy matrix 
into ()A.38p . and taking the normahzed block-trace ()A.28P of both sides, we 
arrive at the 2x2 matrix equation for the elements of the matrix-valued 
Green's function of the Girko-Ginibre ensemble ([3]), 

gzz g^-^ \ _ f z -a^g^^ - ^ ( ^ 

6?^^ g'^ )~\ -a'^g'' z ) ~ |z|2 - a^g^'^g^' \ a^Q'' 

(A.43) 

where in the last step we have explicitly inverted the matrix. Let us first 
look at the negated product of the two off-diagonal equations of (|A.43p . 
which yields an equation for C = Ca{z,z) ()A.33p . 

C = — (A.44) 

{\z\^ + a^Cf 

It has one solution C = 0, corresponding to the outside of V, and called the 
"holomorphic solution," as well as a non-trivial one, describing the inside 
of T), named the "non-holomorphic solution," and given by 

Ca(^,^) = ^(i-^), (A.45) 

where we have picked the positive root out of the two. As mentioned (see the 
discussion below ()A.34p ). setting this quantity to zero yields the equation of 
the borderline of P — since the holomorphic and non-holomorphic solutions 
must coincide on dV, 

\z\ = a, (A.46) 

i.e., the average eigenvalues of the Girko-Ginibre model fill the centered cir- 
cle of radius a. Let us now proceed to the one diagonal equation (stemming 
from the upper left corner) of ()A.43p . 

for l^l < 0", 



r r : (a.47) 



\z\'^ + a^C \ J, for \z\ > a. 

We explicitly see that outside of the eigenvalues' disk, it is trivial, equal 
to the usual Green's function (the "holomorphic solution"), which for the 
Girko-Ginibre ensemble is just 1/z. However, inside the disk, it provides 
the non-holomorphic Green's function ()A.30p of the model, 

Ga{z,z) = ^, (A.48) 

and subsequently (|A.2ip . the mean eigenvalue density (H]). 

The Girko-Ginibre model is therefore solved in the thermodynamic limit 
by means of planar diagrammatics on a "duplicated" level, outlining a pat- 
tern to follow when dealing with more complicated non-Hermitian random 
matrices. 
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Appendix A. 2. 3 Non— Hermitian Planar Diagrammatics for Products of Ran- 
dom Matrices 



As explained in ^2.1.H if one wants to apply the method of planar di- 
agrams to a model being a product of random matrices, the linearization 
(|17p is necessary. In this paragraph, we will prove the main result of this 
procedure, (fT8|l . 

Since the L-th power of ()17p is an L x L block-diagonal matrix (total size 
of course A'tot. x -^tot.)) with the entries being all the cyclic permutations of 
the product ([T]), 



/ A1A2... Al_iAl 
' 





AoAfi . . . Ar Ai 






AlAi ...Ai 



\ 



\ ... AlAi . . . Al-2-^l-i / 

(A.49) 

and because all these permutations have identical non-zero eigenvalues, plus 
a number (distinct for all these matrices) of zero modes — therefore, has 
the same eigenvalues as P, name them Ai, A2, • • • , Aat^ , each one L-fold de- 
generate (so in total LNi), plus (A'tot. — -^-^1) additional zero modes (some 
of the Aa's may^be zero as well). This implies that the content of the A'tot. 
eigenvalues of P is the following: LNi eigenvalues {/%", \/A2, . . . , vX/v^ 
(each one with L-fold degeneracy), where we keep in mind that the L- 
th root is an L-valued function, plus (A'tot. — -^-^1) additional zero modes. 
Thus, we may translate this description into a mathematical expression 
(IA18D . 



a=l 




LNi 

w 



tot. 



5(2) ( 



w, w I 



(A.50) 

Further, we recognize that the piece inside the big brackets in the lower 
line is precisely the mean density of the L-th roots of the eigenvalues of P, 
hence, it will be equal to the distribution of these eigenvalues themselves 
upon the change of variables w = {/z. 



1 ^1 



2) 



a=l 



W 



1 dz dz 
L dw dw 



/9p {W^,W^) 



(A.51) 
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where the factor of 1/L in the Jacobian originates from the L-valuedness 
of the transformation z ^ w. Combining ()A.5ip with ()A.50p . we arrive at 
a relationship between the mean densities of the eigenvalues of P and P, 

p^iw,W) = ^kp(^-^)pp {w^,W^) + ^) 5^'Hw,w). (A.52) 

jVtot. V JVtot. / 

This generalizes formula (16) of [13] to constituent random matrices having 
arbitrary rectangularities. 

Relation (1X521) may be recast (|A.2ip in the language of the non-holomorphic 
Green's function, 



dw ^ dw \Ntot. V ^tot. 

where in the last term on the RHS we have rewritten the complex Dirac's 
delta as 6^'^\w,w) = {l/7T)djjj{l/w), equivalent to ()A.19p upon the regular- 
ization l/w = lim^^QW/{ww + e^). It is of course desirable now to remove 
the derivatives from both sides of (|A.53p . but this could in principle 
produce an additive regular term independent ofW, f{w); however, we will 
shortly justify that f{w) = 0. If so, ()A.53p . without the derivatives, appears 
even simpler in terms of the non-holomorphic M-transform, 

Mp(u;,IZJ) = ^Mp(u;^,^U^), (A.54) 

which is exactly ()18p . 

The idea of proving f{w) = is to show that ()A.54p holds for the 
usual (holomorphic) M-transforms (|A.7p . i.e., outside of the domain T> 
of the average eigenvalues; then, knowing that the non-holomorphic and 
holomorphic M-transforms coincide on the borderline dD, (|A.54p would be 
proven. Hence, expanding the holomorphic M-transforms into the series 
of the moments |A.8p . and noticing that if an integer n is a multiple of L, 
n = Lk, then Tr P" = L Tr P'^ (IXiOjl , while in all the other cases, Tr P" = 
— we obtain 

M^M = V ± J- TV (P"\ = ^ V Tr (P^\ = ^Mp , 

n>l fe>l 

(A.55) 

i.e., precisely like (|A.54p . but in the holomorphic sector. This completes 
the proof of 
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Appendix B 

Free Random Variables in a Nutshell 

In this appendix, we briefly outline some basic facts about the "free 
random variables" (FRV) calculus; given the broadness of the subject, we 
redirect the reader for more details to the specialized literature. The math- 
ematical structure behind FRV, established by Voiculescu et al. and Spe- 
icher [461 147] . is an extension of classical probability theory to the world 
of non-commuting random objects, of which large random matrices are a 
prime example. 

Appendix B.l Freeness and Addition of Random Matrices 
Appendix B.1.1 Independence and the Classical Addition Law 

A fundamental notion in the standard theory of probability is that of 
independence — two random variables Hi and H2 are called "independent" 
if after subtracting their average values, -ff( 2 = ^1,2 — (-^1,2); they fulfil 

{H'lH'^) = 0. (B.l) 

Statistical independence constitutes a basis for numerous important con- 
structions. For instance, the probability distribution of the sum (Hi + H2) 
of independent random numbers depends uniquely on the PDF's of the sum- 
mands, and can be obtained by the following simple procedure: Thanks to 
(jB.ip . any moment of the sum can be expressed through the moments of 
the respective variables using the binomial theorem, 

MH,+H,,n = {{Hi + H^T) = E U {Hlnr') = 

fc=0 ^ ^ 

= E Q (^0 (^2-'^) = E QMn,,^Mn,,n... (B.2) 

A;=0 ^ ^ fc=0 ^ ^ 

In other words, gathering the moments into a generating function, named 
the "characteristic function," 

n>0 \n.>0 / 

where x is a real argument, one discovers that gHi+H2{x) = 9Hi{x)gH2{x)i 
or yet in other words, that the logarithm of the characteristic function, 

ruix) = log gH{x), (B.4) 
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is simply additive when summing independent random numbers, 

rHi+H2{x) = '''Hiix) + rH2ix), for independent ffi, ^^2- (^-5) 
This is the "commutative (classical) addition law." 

Appendix B.1.2 Freeness 

One might wonder whether an analogous structure could be developed 
in the world of random matrices. Let us focus on the Hermitian case first — 
consider L Hermitian N x N random matrices H;, and try to introduce an 
appropriate notion of independence. Let us also take the matrix dimension 

— )• oo, which happens to be a necessary constraint: FRV works only in 
the thermodynamic limit. 

One notices that it is not sufficient to assume statistical independence of 
all the entries of various H/'s for them to be truly called "independent" — 
indeed, one has to additionally require that there is no angular correlation 
between them; such correlations do arise from angular patterns specific to 
the given matrix ensembles. One may remove this angular dependence 
by performing a random similarity transformation of each matrix, with the 
uniform density on the unitary group, U;H;U| — such matrices display then 
all the necessary features, and are referred to as "freely independent," or 
just "free." Non-commutative probability theory equipped with the notion 
of freeness is called "free probability" or "free random variables" (FRV) 
calculus. 

Li order to put this definition in rigorous terms \W\ I47| . one begins 
from making a claim that the appropriate generalization of the classical 
expectation map (. . .) should be r(. . .) = (l/N) Tr(. . .). Now, one says that 
the Hi's are free if their centered counterparts, = H/ — r(H;), obey the 
following analogue of (|B.ip . 

r(H',^H',^...H'J=0, (B.6) 

for all integers n > 2, and all indices /i, • • • i = 1, 2, . . . , L such that 

To get better acquainted with ()B.6p , let us examine some of its simplest 
implications for two free matrices: By writing = + t(H/) and using 
(jB.Gp . one may work out all the correlation functions. The 2-point one is 
then found to be 

r(HiH2) = T(Hi)r(H2), (B.7) 

i.e., precisely like in the commutative case. For the two possible 4-point 
correlation functions, one has 



r (H?H^) = r (H?) r (Hi) , 



(B.8) 
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i.e., again no difference from the scalar sector, but 

r (H1H2H1H2) = r (Hi)' T (Hi) + r (H?) r (H2)' - r (Hi)' r (H2)' , 

(B.9) 

which is very different from the classical case. Hence, mixed moments gener- 
ically do not factorize into moments of separate variables, as in usual prob- 
ability theory — a property which we have seen ()B.2p to be the cornerstone 
of the classical addition law (jB.Sp ; consequently, the latter will require some 
serious amendments to make up for this new situation. 



Appendix B.1.3 The FRV Addition Law for Hermitian Random Matrices 

We have explained ( ^A.l.ip that all the moments of a Hermitian random 
matrix H are conveniently grouped into a generating function such as the 
Green's function Gh(-z) ()A.3p or the M-transform Mh(-2) (|A.7p . In an 
analogous way as in standard probability theory, where the application of 
the logarithm (|B.4p to the characteristic function gnix) ()B.3p leads to the 
object rnix), which is additive (jB.Sp under the addition of independent 
(|B.ip random variables — so in free probability, one has to apply functional 
inversion to the Green's function, 

GH{BH{z)) = BH{Gn{z)) = z, (B.IO) 

finding the so-called "Blue's function," which fulfils the "non-commutative 
(FRV) addition law," 

Bu,+nAz) = BhAz) + BhAz) - I, for free Hi, H2. (B.ll) 

Let us show on an example how this result can be practically used. 
Consider the normalized sum (Hi + H2 + . . . + Hi)/\/L of L free random 
matrices sampled from the Gaussian Unitary Ensembles ()A.9p . with the 
variances cr', (t|, . . . , cr|, respectively. The normalization by is equiv- 

alent to the rescaling of all the variances as erf /L. The Green's function 
of the GUE with an arbitrary variance cr' is given by ()A.16p . It is easy 
to invert it functionally, which produces the corresponding Blue's function 

(iBlnD . 

BHiz)=(T^z + ^. (B.12) 

Hence, the addition law (jB.lip yields the Blue's function of the normalized 
sum to be 

/ N ^ 1 , af + a'^ + ...+a'j 

^(Hi+H2+...+Hi)/VL(^) ='^^ + -' ^'^^'^^ ^ - X • 

(B.13) 
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In other words, it remains a GUE random matrix, yet with the variance 
equal to the arithmetic mean of the constituent variances. 

We have thus observed how one may exploit the FRV addition law ()B.lip 
to make discoveries of non-trivial results concerning sums of freely indepen- 
dent random matrices in an algorithmic, wholly algebraic, and usually quite 
simple way. 

Let us finally remark that there exist several alternative notations com- 
monly used in the literature, equivalent to the Blue's function language 
(jB.lOp . For example, Voiculescu et al. formulated the theory in terms 
of the "/^-transform," RYi{z) = Bii{z) — 1/z, in terms of which the FRV 
addition law (jB.lip appears even simpler, i?Hi+H2(-2) = ^Hi(-2) + Rh2{z) 
(compare (jB.SP ). 

Appendix B.1.4 The FRV Addition Law for Non— Hermitian Random Matrices 

Even though of no use in the current paper, let us mention that the 
FRV addition law (|B.lip can be generalized from the Hermitian to non- 
Hermitian realm with a "quaternion construction" |481 H9] . 

The underlying idea looks as follows: We have elaborated ( ^A.2.ip on 
how the 2x2 matrix-valued Green's function ()A.29p is the proper quan- 
tity to encode all the average spectral information about a non-Hermitian 
matrix model X, 



We have stated that it comprises a non-Hermitian analogue of the Hermitian 
Green's function Gh(-z) (|A.3p . but this is not entirely true — the presence 
of ZP suggests that it is rather a counterpart of Gh(A + ie), e — )• 0"'". On 
the other hand, when one wishes to use the Hermitian FRV addition law 
(jB.lip . one needs the Blue's function, which is the functional inversion of 
the Green's function (jB.lOp . and hence, the functional form of the latter is 
necessary for any complex z, and not only for z = \ + ie. 

In a pursuit of a non-Hermitian generalization of the Blue's function, one 
is therefore compelled to replace ZP in ()B.14p — sufficient for producing the 
average spectral density, yet incapable of engineering functional inversion 
— by an arbitrary quaternion, 



where a, b are complex numbers. (Actually, one might restrain b to the real 
axis, i.e., choose a 3-dimensional subspace in the quaternion space. This, 
however, does not result in simplifying any equation, so we will abandon it. 




(B.14) 




(B.15) 
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The point is that b must no longer be close to zero, but cover the entire real 
line/complex plane.) In other words, the matrix-valued Green's function 
(|B.14p is accordingly promoted to a quaternion function of a quaternion 
variable, named the "quaternion Green's function," 

gx(Q) = ^lini^^bTr((Q®ljv-XD)"y (B.16) 

After this modification, functional inversion is finally allowed in the 
quaternion space, and one defines the "quaternion Blue's function," 

(ex(Q)) = (ax(Q)) = Q. (B.17) 

Remarkably H9] , it may be proven that it satisfies the "quaternion FRV 
addition law," which parallels the Hermitian version (IRTT]) . but at the 
quaternion level, 

fixi+X2(Q) = eXi(Q) + ex2(Q)-^, for free Xi, X2. (B.18) 

This result constitutes a robust tool, greatly simplifying handling of sums 
of free non-Hermitian random matrices (see for instance |5m to view it 
at work). 

Appendix B. 2 Multiplication of Random Matrices 

Another important problem, after the addition of free random matri- 
ces, is their multiplication. For scalar random variables, it boils down 
to the addition task through the exponential change of variables, since 
exp Hi exp H2 = exp{Hi + H2). Unfortunately, for non-commuting objects, 
generically exp Hi exp H2 7^ exp(Hi + H2), and one faces a challenge to de- 
vise a different approach. 

The FRV calculus provides such a straightforward prescription of mul- 
tiplying free random matrices. They must be either Hermitian, but with 
an additional assumption that their product is Hermitian, too, or unitary, 
in which case the product automatically remains unitary. The procedure 
for such matrices looks as follows: First, invert functionally the pertinent 
M-transforms (|A.7p . 

Mh (A^h(^)) = Nn {Mn{z)) = z, (B.19) 

obtaining the so-named "A^-transforms." This object is proven to satisfy 
the "non-commutative (FRV) multiplication law," 

Nu.nA^) = ^-NnAz)Nn,iz), for free Hi, Hg. (B.20) 
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We exploit this formula extensively in section [3l 

Finally, let us mention that our notation is distinct from the one most 
often used in the context of multiplying free random variables; one may 
alternatively define the "S-transform," S'h(-z) = {z + 1) / {zN-f:i{z)) , which 
is simply multiplicative, S-h.^-h_^{z) = S-h.^{z)S-h.^[z). 
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